if(FALSE){
  #setwd("C:\\Research\\csiro\\Species13\\Ling\\2013 Runs\\Spatial1\\")
  #setwd('c:/SS/Thornyheads/runs/SST_base_FINAL/')

  setwd("C:/ss/modeltesting/Version_3_24s_July24/Simple")
  CompData <- read.table("CompReport.SSO",fill=T,skip=21,head=T,as.is=T)
  AgeData <- CompData[CompData$Kind == "AGE",]
  AgeData$Yr <- as.double(AgeData$Yr)
  AgeData$Gender <- as.double(AgeData$Gender)
  LengthData <- CompData[CompData$Kind == "LEN",]
  LengthData$Yr <- as.double(LengthData$Yr)
}

LengthCase <- function()
{
 Fleets <- sort(unique(LengthData$Fleet))
 print(Fleets)
 par(mfrow=c(2,3))
 for (Ifleet in Fleets)
  {
   TheD <- LengthData[LengthData$Fleet== Ifleet,] 
   Years <- sort(unique(TheD$Yr))
   YrsV <- rep(0,length(Years))
   ObsV <- rep(0,length(Years))
   ExpE <- rep(0,length(Years))
   VarE <- rep(0,length(Years))
   Ipnt <- 0
   for (Yr in Years)
    {
     Ipnt <- Ipnt + 1
     TheD2 <- TheD[TheD$Yr==Yr & !is.na(TheD$Bin),] 
     EffN <- sum(TheD2$Exp*(1-TheD2$Exp))/sum((TheD2$Obs-TheD2$Exp)^2)
     MidLen <- TheD2$Bin + 2.5
     ObsV[Ipnt] <- sum(TheD2$Obs*MidLen)
     ExpE[Ipnt] <- sum(TheD2$Exp*MidLen)
     N <- TheD2$N[1]
     YrsV[Ipnt] <- Yr
     VarE[Ipnt] <- sqrt((sum(TheD2$Exp*MidLen^2) - ExpE[Ipnt]^2)/N)
    }
   Resid <- (ObsV-ExpE)/VarE
   Adjust = 1/var(Resid)
   #print(ObsV)  
   #print(ExpE)
   #print(VarE)
   #print(YrsV)
   PredV1Low <- ObsV-VarE*1.96
   PredV1Hi <- ObsV+VarE*1.96
   VarE2 <- VarE/sqrt(Adjust)
   PredV2Low <- ObsV-VarE2*1.96
   PredV2Hi <- ObsV+VarE2*1.96
  
   ymax <- max(ObsV,PredV1Hi,PredV2Hi)*1.1
   ymin <- min(ObsV,PredV1Low,PredV2Low)/1.1
   plot(0,0,xlab="Year",ylab="Mean Length",xlim=c(1977,2012),ylim=c(ymin,ymax))
   points(YrsV,ObsV,pch=16)
   lines(YrsV,ExpE,lty=1,lwd=2)
   for (II in 1:Ipnt)
    {
     lines(c(YrsV[II],YrsV[II]),c(PredV1Low[II],PredV1Hi[II]),lty=1,lwd=2)
     lines(c(YrsV[II],YrsV[II]),c(PredV2Low[II],PredV2Hi[II]),lty=1,lwd=1)
    }  
   cat("Length adjustment for fleet ",Ifleet," = ",round(Adjust,3),"\n")
   title(paste("Fleet =",Ifleet))
  }  
}

AgeCase <- function(Diag=F)
{
  Fleets <- sort(unique(AgeData$Fleet))
  par(mfrow=c(2,3))
  for (Ifleet in Fleets)
  {
    TheD <- AgeData[AgeData$Fleet==Ifleet,] 
    Years <- sort(unique(TheD$Yr))
    Resid <- NULL
    for (Yr in Years)
     for (Gender in 0:2)  
     {
      Use <-TheD$Yr==Yr & TheD$Pick_gender==Gender & !is.na(TheD$Bin)
      if (Diag==T & sum(Use) == 0) cat(Yr,Gender,sum(Use),"\n")
      if (sum(Use) > 0)
       {
        TheD2 <- TheD[Use,]
        LowBin <- unique(TheD2$Lbin_lo)
        #cat(Yr,Gender,sum(Use)/21,LowBin,"\n")
        for (bins in LowBin)
         {
          TheD3 <- TheD2[TheD2$Lbin_lo == bins,]
          ObsVA <- sum(TheD3$Obs*TheD3$Bin)
          ExpEA <- sum(TheD3$Exp*TheD3$Bin)
          N <- TheD3$N[1]
          if (N == 0) { print("oops"); AAA}
          VarEA <- sqrt((sum(TheD3$Exp*TheD3$Bin^2) - ExpEA^2)/N)
          if (Diag == T) cat(ObsVA,ExpEA,VarEA,"\n")
          Resid <- c(Resid,(ObsVA-ExpEA)/VarEA)
         }  
      }    
     }  
    #print(Resid)
    hist(Resid,main="")
    title(paste("Fleet=",Ifleet))
    Adjust = 1/var(Resid)
    cat("Age adjustment for fleet ",Ifleet," = ",round(Adjust,3),"\n")
  }  
}

if(FALSE){
  LengthCase()
  AgeCase()
}
